home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick_docs folder / yorick_docs / math.doc < prev    next >
Text File  |  1996-02-28  |  19KB  |  440 lines
  1.  
  2.  
  3.  
  4.                           Yorick Documentation
  5.                 for functions, variables, and structures
  6.          defined in files fft.i and matrix.i
  7.                    Printed: Sat Mar 18 22:52:20 1995
  8.  
  9.    Contents:
  10.  
  11.                                    LUrcond
  12.                                    LUsolve
  13.                                    QRsolve
  14.                                    SVdec
  15.                                    SVsolve
  16.                                    TDsolve
  17.                                    _dgecox
  18.                                    _dgelss
  19.                                    _dgelx
  20.                                    _dgesv
  21.                                    _dgesvx
  22.                                    _dgetrf
  23.                                    _dgtsv
  24.                                    fft
  25.                                    fft_braw
  26.                                    fft_fraw
  27.                                    fft_init
  28.                                    fft_inplace
  29.                                    fft_setup
  30.                                    roll
  31.                                    unit
  32.  
  33.                                                     FROM LUrcond TO LUsolve
  34.  
  35.                                                                     LUrcond
  36. /* DOCUMENT LUrcond(a)
  37.          or LUrcond(a, one_norm=1)
  38.      returns the reciprocal condition number of the N-by-N matrix A.
  39.      If the ONE_NORM argument is non-nil and non-zero, the 1-norm
  40.      condition number is returned, otherwise the infinity-norm condition
  41.      number is returned.
  42.  
  43.      The condition number is the ratio of the largest to the smallest
  44.      singular value, max(singular_values)*max(1/singular_values) (or
  45.      sum(abs(singular_values)*sum(abs(1/singular_values)) if ONE_NORM
  46.      is selected?).  If the reciprocal condition number is near zero
  47.      then A is numerically singular; specifically, if
  48.           1.0+LUrcond(a) == 1.0
  49.      then A is numerically singular.
  50.  
  51.    SEE ALSO: LUsolve
  52.  */
  53.  
  54.                                                                     LUsolve
  55. /* DOCUMENT LUsolve(a, b)
  56.          or LUsolve(a, b, which=which)
  57.      or a_inverse= LUsolve(a)
  58.  
  59.      returns the solution x of the matrix equation:
  60.         A(,+)*x(+) = B
  61.      If A is an n-by-n matrix then B must have length n, and the returned
  62.      x will also have length n.
  63.  
  64.      B may have additional dimensions, in which case the returned x
  65.      will have the same additional dimensions.  The WHICH dimension of B,
  66.      and of the returned x is the one of length n which participates
  67.      in the matrix solve.  By default, WHICH=1, so that the equations
  68.      being solved are:
  69.         A(,+)*x(+,..) = B
  70.      Non-positive WHICH counts from the final dimension (as for the
  71.      sort and transpose functions), so that WHICH=0 solves:
  72.         x(..,+)*A(,+) = B
  73.  
  74.      If the B argument is omitted, the inverse of A is returned:
  75.      A(,+)*x(+,) and A(,+)*x(,+) will be unit matrices.
  76.  
  77.      LUsolve works by LU decomposition using Gaussian elimination with
  78.      pivoting.  It is the fastest way to solve square matrices.  QRsolve
  79.      handles non-square matrices, as does SVsolve.  SVsolve is slowest,
  80.      but can deal with highly singular matrices sensibly.
  81.  
  82.    SEE ALSO: QRsolve, TDsolve, SVsolve, SVdec, LUrcond
  83.  */
  84.  
  85.                                                       FROM QRsolve TO SVdec
  86.  
  87.                                                                     QRsolve
  88. /* DOCUMENT QRsolve(a, b)
  89.          or QRsolve(a, b, which=which)
  90.  
  91.      returns the solution x (in a least squares or least norm sense
  92.      described below) of the matrix equation:
  93.         A(,+)*x(+) = B
  94.      If A is an m-by-n matrix (i.e.- m equations in n unknowns), then B
  95.      must have length m, and the returned x will have length n.
  96.  
  97.      If n<m, the system is overdetermined -- no solutions are possible
  98.              -- the returned x minimizes sqrt(sum((A(,+)*x(+) - B)^2))
  99.      If n>m, the system is underdetermined -- many solutions are possible
  100.              -- the returned x has minimum L2 norm among all solutions
  101.  
  102.      B may have additional dimensions, in which case the returned x
  103.      will have the same additional dimensions also have those dimensions.
  104.      The WHICH dimension of B and the returned x is the one of length m
  105.      or n which participates in the matrix solve.  By default, WHICH=1,
  106.      so that the equations being solved are:
  107.         A(,+)*x(+,..) = B
  108.      Non-positive WHICH counts from the final dimension (as for the
  109.      sort and transpose functions), so that WHICH=0 solves:
  110.         A(,+)*x(..,+) = B
  111.  
  112.      QRsolve works by QR factorization if n<m, or LQ factorization if n>m.
  113.      QRsolve is slower than LUsolve.  Its main attraction is that it can
  114.      handle overdetermined or underdetermined systems of equations
  115.      (nonsquare matrices).  QRsolve may fail for singular systems; try
  116.      SVsolve in this case.
  117.  
  118.    SEE ALSO: LUsolve, TDsolve, SVsolve, SVdec
  119.  */
  120.  
  121.                                                                       SVdec
  122. /* DOCUMENT s= SVdec(a, u, vt)
  123.          or s= SVdec(a, u, vt, full=1)
  124.  
  125.      performs the singular value decomposition of the m-by-n matrix A:
  126.         A = (U(,+) * SIGMA(+,))(,+) * VT(,+)
  127.      where U is an m-by-m orthogonal matrix, VT is an n-by-n orthogonal
  128.      matrix, and SIGMA is an m-by-n matrix which is zero except for its
  129.      min(m,n) diagonal elements.  These diagonal elements are the return
  130.      value of the function, S.  The returned S is always arranged in
  131.      order of descending absolute value.  U(,1:min(m,n)) are the left
  132.      singular vectors corresponding to the min(m,n) elements of S;
  133.      VT(1:min(m,n),) are the right singular vectors.  (The original A
  134.      matrix maps a right singular vector onto the corresponding left
  135.      singular vector, stretched by a factor of the singular value.)
  136.  
  137.      Note that U and VT are strictly outputs; if you don't need them,
  138.      they need not be present in the calling sequence.
  139.  
  140.      By default, U will be an m-by-min(m,n) matrix, and V will be
  141.      a min(m,n)-by-n matrix (i.e.- only the singular vextors are returned,
  142.  
  143.                                                       FROM SVdec TO SVsolve
  144.  
  145.      not the full orthogonal matrices).  Set the FULL keyword to a
  146.      non-zero value to get the full m-by-m and n-by-n matrices.
  147.  
  148.      On rare occasions, the routine may fail; if it does, the
  149.      first SVinfo values of the returned S are incorrect.  Hence,
  150.      the external variable SVinfo will be 0 after a successful call
  151.      to SVdec.  If SVinfo>0, then external SVe contains the superdiagonal
  152.      elements of the bidiagonal matrix whose diagonal is the returned
  153.      S, and that bidiagonal matrix is equal to (U(+,)*A(+,))(,+) * V(+,).
  154.  
  155.      Numerical Recipes (Press, et. al. Cambridge University Press 1988)
  156.      has a good discussion of how to use the SVD -- see section 2.9.
  157.  
  158.    SEE ALSO: SVsolve, LUsolve, QRsolve, TDsolve
  159.  */
  160.  
  161.                                                                     SVsolve
  162. /* DOCUMENT SVsolve(a, b)
  163.          or SVsolve(a, b, rcond)
  164.          or SVsolve(a, b, rcond, which=which)
  165.  
  166.      returns the solution x (in a least squares sense described below) of
  167.      the matrix equation:
  168.         A(,+)*x(+) = B
  169.      If A is an m-by-n matrix (i.e.- m equations in n unknowns), then B
  170.      must have length m, and the returned x will have length n.
  171.  
  172.      If n<m, the system is overdetermined -- no solutions are possible
  173.              -- the returned x minimizes sqrt(sum((A(,+)*x(+) - B)^2))
  174.      If n>m, the system is underdetermined -- many solutions are possible
  175.              -- the returned x has minimum L2 norm among all solutions
  176.  
  177.      SVsolve works by singular value decomposition, therefore it is
  178.      immune to failure due to singularity of the A matrix.  The optional
  179.      RCOND argument defaults to 1.0e-9; singular values less than RCOND
  180.      times the largest singular value (absolute value) will be set to 0.0.
  181.      If RCOND<=0.0, machine precision is used.  The effective rank of the
  182.      matrix is returned as the external variable SVrank.
  183.  
  184.      You can examine the details of the SVD by calling the SVdec routine,
  185.      which returns the singular vectors as well as the singular values.
  186.      Numerical Recipes (Press, et. al. Cambridge University Press 1988)
  187.      has a good discussion of how to use the SVD -- see section 2.9.
  188.  
  189.      B may have additional dimensions, in which case the returned x
  190.      will have the same additional dimensions.  The WHICH argument
  191.      (default 1) controls which dimension of B takes part in the matrix
  192.      solve.  See QRsolve or LUsolve for a complete discussion.
  193.  
  194.    SEE ALSO: SVdec, LUsolve, QRsolve, TDsolve
  195.  */
  196.  
  197.                                                     FROM TDsolve TO _dgesvx
  198.  
  199.                                                                     TDsolve
  200. /* DOCUMENT TDsolve(c, d, e, b)
  201.          or TDsolve(c, d, e, b, which=which)
  202.  
  203.      returns the solution to the tridiagonal system:
  204.         D(1)*x(1)       + E(1)*x(2)                       = B(1)
  205.     C(1:-1)*x(1:-2) + D(2:-1)*x(2:-1) + E(2:0)*x(3:0) = B(2:-1)
  206.                       C(0)*x(-1)      + D(0)*x(0)     = B(0)
  207.      (i.e.- C is the subdiagonal, D the diagonal, and E the superdiagonal;
  208.      C and E have one fewer element than D, which is the same length as
  209.      both B and x)
  210.  
  211.      B may have additional dimensions, in which case the returned x
  212.      will have the same additional dimensions.  The WHICH dimension of B,
  213.      and of the returned x is the one of length n which participates
  214.      in the matrix solve.  By default, WHICH=1, so that the equations
  215.      being solved involve B(,..) and x(+,..).
  216.      Non-positive WHICH counts from the final dimension (as for the
  217.      sort and transpose functions), so that WHICH=0 involves B(..,)
  218.      and x(..,+).
  219.  
  220.      The C, D, and E arguments may be either scalars or vectors; they
  221.      will be broadcast as appropriate.
  222.  
  223.   SEE ALSO: LUsolve, QRsolve, SVsolve, SVdec
  224.  */
  225.  
  226.                                                                     _dgecox
  227. /* DOCUMENT _dgecox
  228.      LAPACK dgecon routine, except norm argument not a string.
  229.  */
  230.  
  231.                                                                     _dgelss
  232. /* DOCUMENT _dgelss
  233.      LAPACK dgelss routine.
  234.  */
  235.  
  236.                                                                      _dgelx
  237. /* DOCUMENT _dgelx
  238.      LAPACK dgels routine, except trans argument not a string.
  239.  */
  240.  
  241.                                                                      _dgesv
  242. /* DOCUMENT _dgesv
  243.      LAPACK dgesv routine.
  244.  */
  245.  
  246.                                                                     _dgesvx
  247. /* DOCUMENT _dgesvx
  248.      LAPACK dgesvd routine, except jobu and jobvt are not strings.
  249.  */
  250.  
  251.                                                         FROM _dgetrf TO fft
  252.  
  253.                                                                     _dgetrf
  254. /* DOCUMENT _dgetrf
  255.      LAPACK dgetrf routine.  Performs LU factorization.
  256.  */
  257.  
  258.                                                                      _dgtsv
  259. /* DOCUMENT _dgtsv
  260.      LAPACK dgtsv routine.
  261.  */
  262.  
  263.                                                                         fft
  264. /* DOCUMENT fft(x, direction)
  265.             fft(x, ljdir, rjdir)
  266.          or fft(x, ljdir, rjdir, setup=workspace)
  267.      returns the complex Fast Fourier Transform of array X.
  268.      The DIRECTION determines which direction the transform is in --
  269.      e.g.- from time to frequency or vice-versa -- as follows:
  270.  
  271.      DIRECTION    meaning
  272.      ---------    -------
  273.          1        "forward" transform (coefficients of exp(+i * 2*pi*kl/N))
  274.               on every dimension of X
  275.         -1        "backward" transform (coefficients of exp(-i * 2*pi*kl/N))
  276.               on every dimension of X
  277.      [1,-1,1]     forward transform on first and third dimensions of X,
  278.                   backward transform on second dimension of X (any other
  279.           dimensions remain untransformed)
  280.      [-1,0,0,1]   backward transform on first dimension of X, forward
  281.                   transform on fourth dimension of X
  282.     etc.
  283.  
  284.      The third positional argument, if present, allows the direction
  285.      of dimensions of X to be specified relative to the final dimension
  286.      of X, instead of relative to the first dimension of X.  In this
  287.      case, both LJDIR and RJDIR must be vectors of integers -- the
  288.      scalar form is illegal:
  289.  
  290.         LJDIR    RJDIR      meaning
  291.         -----    -----      -------
  292.     []        [1]       forward transform last dimension of X
  293.     [1]        []       forward transform first dimension of X
  294.     []        [-1,-1]   backward transform last two dimensions of X,
  295.                         leaving any other dimensions untransformed
  296.      [-1,0,0,1]    []       backward transform on first dimension of X,
  297.                             forward transform on fourth dimension of X
  298.         []      [-1,0,0,1]  backward transform on 4th to last dimension of X,
  299.                             forward transform on last dimension of X
  300.     etc.
  301.  
  302.      Note that the final element of RJDIR corresponds to the last dimension
  303.      of X, while the initial element of LJDIR corresponds to the first
  304.      dimension of X.
  305.  
  306.      The explicit meaning of "forward" transform -- the coefficients of
  307.      exp(+i * 2*pi*kl/N) -- is:
  308.  
  309.                                                     FROM fft TO fft_inplace
  310.  
  311.  
  312.      result for j=1,...,n
  313.  
  314.                 result(j)=the sum from k=1,...,n of
  315.  
  316.                       x(k)*exp(-i*(j-1)*(k-1)*2*pi/n)
  317.  
  318.                             where i=sqrt(-1)
  319.  
  320.      Note that the result is unnormalized.  Applying the "backward"
  321.      transform to the result of a "forward" transform returns N times
  322.      the original vector of length N.  Equivalently, applying either
  323.      the "forward" or "backward" transform four times in succession
  324.      yields N^2 times the original vector of length N.
  325.  
  326.      Performing the transform requires some WORKSPACE, which can be
  327.      set up beforehand by calling fft_setup, if fft is to be called
  328.      more than once with arrays X of the same shape.  If no setup
  329.      keyword argument is supplied, the workspace allocation and setup
  330.      must be repeated for each call.
  331.  
  332.    SEE ALSO: roll, fft_setup, fft_inplace
  333.  */
  334.  
  335.                                                                    fft_braw
  336. /* DOCUMENT fft_braw, n, c, wsave
  337.      Swarztrauber's cfftb.  You can use this to avoid the additional
  338.      2*N storage incurred by fft_setup.
  339.  */
  340.  
  341.                                                                    fft_fraw
  342. /* DOCUMENT fft_fraw, n, c, wsave
  343.      Swarztrauber's cfftf.  You can use this to avoid the additional
  344.      2*N storage incurred by fft_setup.
  345.  */
  346.  
  347.                                                                    fft_init
  348. /* DOCUMENT fft_init, n, wsave
  349.      Swarztrauber's cffti.  This actually requires wsave=array(0.0, 4*n+15),
  350.      instead of the 6*n+15 doubles of storage used by fft_raw to handle the
  351.      possibility of multidimensional arrays.  If the storage matters, you
  352.      can call cfftf and/or cfftb as the Yorick functions fft_fraw and/or
  353.      fft_braw.
  354.  */
  355.  
  356.                                                                 fft_inplace
  357. /* DOCUMENT fft_inplace, x, direction
  358.          or fft_inplace, x, ljdir, rjdir
  359.      or fft_inplace, x, ljdir, rjdir, setup=workspace
  360.      is the same as the fft function, except that the transform is
  361.      performed "in_place" on the array X, which must be of type complex.
  362.    SEE ALSO: fft, fft_setup
  363.  */
  364.  
  365.                                                      FROM fft_setup TO roll
  366.  
  367.                                                                   fft_setup
  368. /* DOCUMENT workspace= fft_setup(dimsof(x))
  369.          or workspace= fft_setup(dimsof(x), direction)
  370.          or workspace= fft_setup(dimsof(x), ljdir, rjdir)
  371.      allocates and sets up the workspace for a subsequent call to
  372.             fft(X, DIRECTION, setup=WORKSPACE)
  373.      or
  374.             fft(X, LJDIR, RJDIR, setup=WORKSPACE)
  375.      The DIRECTION or LJDIR, RJDIR arguments compute WORKSPACE only for
  376.      the dimensions which will actually be transformed.  If only the
  377.      dimsof(x) argument is supplied, then WORKSPACE will be enough to
  378.      transform any or all dimensions of X.  With DIRECTION or LJDIR, RJDIR
  379.      supplied, WORKSPACE will only be enough to compute the dimensions
  380.      which are actually to be transformed.  The WORKSPACE does not
  381.      depend on the sign of any element in the DIRECTION (or LJDIR, RJDIR),
  382.      so you can use the same WORKSPACE for both "forward" and "backward"
  383.      transforms.
  384.  
  385.      Furthermore, as long as the length of any dimensions of the array
  386.      X to be transformed are present in WORKSPACE, it may be used in
  387.      a call to fft with the array.  Thus, if X were a 25-by-64 array,
  388.      and Y were a 64-vector, the following sequence is legal:
  389.           ws= fft_setup(dimsof(x));
  390.       xf= fft(x, 1, setup=ws);
  391.       yf= fft(y, -1, setup=ws);
  392.  
  393.      The WORKSPACE required for a dimension of length N is 6*N+15 doubles.
  394.  
  395.    SEE ALSO: fft, dimsof, fft_inplace
  396.  */
  397.  
  398.                                                                        roll
  399. /* DOCUMENT roll(x, ljoff, rjoff)
  400.          or roll, x, ljoff, rjoff
  401.      or roll(x)
  402.      or roll, x
  403.  
  404.      "rolls" selected dimensions of the array X.  The roll offsets
  405.      LJOFF and RJOFF (both optional) work in the same fashion as the
  406.      LJDIR and RJDIR arguments to the fft function:
  407.  
  408.         A scalar LJDIR (and nil RJDIR) rolls all dimensions of X by
  409.         the specified offset.
  410.     Otherwise, the elements of the LJDIR vector [ljoff1, ljoff2, ...]
  411.     are used as the roll offsets for the first, second, etc.
  412.     dimensions of X.
  413.     Similarly, the elements of the RJDIR vector [..., rjoff1, rjoff0]
  414.     are matched to the final dimensions of X, so the next to last
  415.     dimension is rolled by rjoff1 and the last dimension by rjoff0.
  416.  
  417.     As a special case (mostly for use with the fft function), if
  418.     both LJDIR and RJDIR are nil, every dimension is rolled by
  419.     half of its length.  Thus,
  420.        roll(x)
  421.     it equivalent to
  422.  
  423.                                                           FROM roll TO unit
  424.  
  425.        roll(x, dimsof(x)(2:0)/2)
  426.  
  427.      The result of the roll function is complex if X is complex, otherwise
  428.      double (i.e.- any other array type is promoted to type double).  If
  429.      roll is invoked as a subroutine, the operation is performed in place.
  430.  
  431.    SEE ALSO: fft
  432.  */
  433.  
  434.                                                                        unit
  435. /* DOCUMENT unit(n)
  436.          or unit(n, m)
  437.      returns n-by-n (or n-by-m) unit matrix, i.e.- matrix with diagonal
  438.      elements all 1.0, off diagonal elements 0.0
  439.  */
  440.